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Abstract. We provide an exact expression for the statistics of the fluxes of Markov 
jump processes at all times, improving on asymptotic results from large deviation 
theory. The main ingredient is a generalization of the BEST theorem in enumeratoric 
graph theory to Eulerian tours with open ends. In the long-time limit we reobtain 
Sanov’s theorem for Markov processes, which expresses the exponential suppression of 
fluctuations in terms of relative entropy. The finite-time power-law term, increasingly 
important with the system size, is a spanning-tree determinant that, by introducing 
Grassmann variables, can be absorbed into the effective Lagrangian of a Fermionic 
ghost field on a metric space, coupled to a gauge potential. With reference to concepts 
in nonequilibrium stochastic thermodynamics, the metric is related to the dynamical 
activity that measures net communication between states, and the connection is made 
to a previous gauge theory for diffusion processes. 
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1. Introduction 

One of the lesser-known intuitions of Einstein, which he applied the critical opalescence 
of mixtures, was the idea of “inverting Boltzmann’s formula” to obtain an expression 
for the probability of a fluctuating variable in terms of entropy jlj. A similar mechanism 
underlies the important result by Sanov [2] in what came to be known as the theory 
of large deviations (LDT): relative entropy describes the asymptotic behavior of the 
probability of an increasing number of independent and identically distributed variables. 
The theory has since grown to embrace Markov processes and beyond. 

Recent years have witnessed a burst of interest in the characterization of the 
thermodynamics of nonequilibrium systems by means of classical and quantum Markov 
processes [3|. The celebrated Fluctuation Theorem |U |5j f6, IQj, which expresses 

a symmetry of the large fluctuations of thermodynamic observables, paved the way for 
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a thorough application of the full spectrum of methods and results from LDT m- 
Here, “large” often indicates the long time limit, regulating the exposure of the system 
to external agents, whereas in equilibrium statistical mechanics one usually refers to 
size. This revival brought about computational mm and theoretical mm m m 
techniques for a more efficient computation of rate functions and cumulants of time- 
extensive observables, a better understanding of the analogies and differences between 
equilibrium and nonequilibrium systems da EH EH EU, of phase transitions along most 
probable paths [221123]- of the hydrodynamic limit fMj and of exotic properties of time¬ 
intensive observables (e.g. the efficiency) |25l [26] . 

Large deviation principles meliorate the central limit theorem. Still there are further 
margins for improvement in considering finite-time effects. However, unfortunately, 
obtaining rate functions in LDT is often prohibitive, but for very simple or special 
systems EH]- The reason for this is that macroscopic observables of interest are 
aggregate quantities to whom contribute many, more fundamental, correlated processes. 
General expressions are only available when one refines the description to all empirical 
fluxes and densities populating the system [29]. 

Indeed, as we will prove in this paper, for Markov processes on a discrete state 
space, this level of description (sometimes called “level 2.5”) is so complete that it 
actually allows for the derivation of the full joint probability of the fluxes and of the 
empirical measure. The finite-time corrections are accounted for by a combinatorial 
factor counting Eulerian tours or, equivalently, spanning trees on a suitably defined 
digraph. Employing the same tools as in a recently proposed Fermionic field theory for 
trees and forests [30l [3Tj . we express the new finite-time term as an effective Lagrangian 
of a ghost held coupled to a gauge potential, on a metric. The metric depends on 
the dynamical activity of the system [29], and a connection is made to a previously 
proposed gauge theory of nonequilibrium thermodynamics by the Author [32] . Finally, 
we re-derive large deviation rate functions by a rather different procedure than is usually 
employed. This direct evaluation allows to appreciate a subtlety related to the large 
deviations of the statistics of the total number of jumps, which is a stochastic variable 
affecting the set of measurable events along a trajectory. 

To facilitate the reading, we will first provide all results in Sec. [2] in the simpler 
setting of a restless, discrete-time Markov chain, reobtaining Sanov’s theorem for 
Markov chains and discussing finite-time corrections to small fluctuations. In Sec. [3] 
we discuss the ghost Lagrangian, and in Sec. [4] we generalize to Markov jump processes 
in continuous time. When considering large deviations for the joint fluxes and empirical 
densities, in § 14.41 we discuss the subtle role played by the total number of jumps. We 
draw conclusions in Sec. [5] 
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2. The simplest case: restless Markov chains 

2.1. Setup 

Consider a graph G = (X,E) consisting of a finite number of vertices x,y,... e X 
connected by edges xy e E. We suppose the graph has no loops (edges that connect 
a state to itself) nor multiple edges. The graph is undirected (arrows do not have 
a preferential direction), but eventually an arbitrary orientation can be established, 
which will be useful in Sec. [3l 

On such state space, a Markovian walker performs a path x = (x 0 ,... ,xn), 
departing from xq and at intervals of unit time performing transitions to neighboring 
vertices Xj. If edge xy exists, then both the probability n xy of jumping from y to x and 
its reverse n yx are nonvanishing, while we assume null probability of resting at states, 
tt xx = 0, Vx. The probability of a path x is given by 

N -1 

P[x] = n TT Xi+lXi . (1) 

2=0 

Introducing the transition matrix II = (n xy ) x . y , the probability pflx) of being at x at 
time i satisfies p i+1 = lip*, given po(x) = 5 XjX0 . Under the above assumptions there exists 
a unique invariant probability p that satisfies lip = p = pjv-*•+«>■ 

The objects of crucial interest in this paper are the fluxes, stochastic variables that 
count the number of transitions between states along a path: 

N -1 

■- X! ^Xi+i,x^Xi,y (2) 

2=0 

Flux <p xy increases by +1 every time transition x <- y occurs. We collectively denote 
the fluxes <p = (p xy ) Xi y. Notice that Y, x ,y wM = -N. It is often physically meaningful 
to focus on the symmetric and antisymmetric parts of the fluxes, called respectively 
activities and currents: 

Txy ~ 2 (dPxy ± Wyx) ■ ( 3 ) 

Physical observables typically have the form 'Z x , y obs xy (p X y[x\. In particular in the 
context of nonequilibrium statistical mechanics, thermodynamic observables like the 
entropy production, the heat fluxes etc. are antisymmetric, obs^ = -obs ya: , hence their 
statistics only depends on the currents. However, while such antisymmetric observables 
suffice to capture the physics near equilibrium, e.g. via the fluctuation-dissipation 
relation and the minimum entropy production principle [33], it is now understood that 
out of equilibrium also symmetric observables related to dynamical effects enter into 
play in a crucial manner. For example, they are necessary to extend the fluctuation- 
dissipation relation out of nonequilibrium steady states ( 34] . 

Importantly, any such kind of observable only depends on the fluxes and not on the 
precise trajectory that generates that configuration of fluxes. Then a preliminary step 


BEST statistics of Markovian fluxes: a tale of Eulerian tours and Fermionic ghosts 4 
is to obtain the probability of the fluxes 

Plfl'EWfl'l ( 4 ) 

X 

= J V(...) (5) 

where <5 is a product of Kronecker deltas, one per flux, and the second expression is a 
sort of yet-to-determine representation in terms of an effective Lagrangian. Asymptotic 
expressions for N -*■ +oo are well-known. The task of this paper is to exactly evaluate 
PW\ at all times and to make sense of the second expression. The problem of 
further marginalizing to evaluate the statistics of actual observables is intractable in 
full generality. 


2.2. BEST statistics of the fluxes 

Given Eq. (J2J) , we can recast the path probability Eq. (pQ) as 


P[x\ = n 


7T 


ifxy 
xy ■ 


x,y 


We can then marginalize for the fluxes to obtain 


P[<P] = £ \.^P\ ex P ( X! Txy In 


t r 


xy I i 


\x,y 


where 


A<P\ = E S 




( 6 ) 


(7) 


( 8 ) 


is the number of distinguishable paths that perform tp xy times transition x •«- y, for all 
transitions. 

The key step in our treatment is the determination of this number, basing on a 
generalization of the so-called BESdljj theorem in enumeratoric graph theory. Let us 
describe the general procedure; in |Appendix C we provide an example of what follows. 
Corresponding to a set of fluxes ip = (p xy ) x , y we can draw a digraph (X,ip) with ip xy 
arrows y -*■ x. In this case, the arrows are considered to be directed. The influx and the 
outflux (in graph-theoretic jargon, in-degree and out-degree ) at vertex x are respectively 
defined as the number of transitions towards and out of x, 

Tx[ x ] : ='Z ( Txy[x], T x ut i x ] (9) 

y y 

Since a Markovian path can be drawn along the edges of the graph by “never lifting the 
pencil”, then the influxes and the outfluxes of all vertices but the initial and the final 
ones balance each other, and in general we have 


T x [x] - p x ut [x] = 5 ; 


xn,x 


j Xq,X‘ 


( 10 ) 


When xq - xN the digraph is said to be balanced. Historically the first theorem in 
graph theory, motivated by the problem whether it was possible to traverse the seven 


f After de Bruijn, van Aardenne-Ehrenfest, Smith and Tutte. 
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Konisberg bridges exactly once to arrive at the place of departure, was Euler’s proof 
that a digraph is balanced if and only if there exists a cyclic path, called Eulerian tour , 
that traverses each arrow exactly once. Importantly for what follows, notice that in an 
Eulerian tour on a balanced digraph each arrow is accounted distinctly, while the paths 
we are considering are blind to the arrow’s identity. The BEST theorem in enumeratoric 
graph theory counts the number of different Eulerian tours on a balanced digraph (see 

m, §5.6). 


Generalizing, we will call a digraph whose in- and out- degrees obey Eq. (fTOl) a quasi- 


balanced digraph. Vertex x 0 is called the source and vertex the sink. In Appendix A 


it is proven that a digraph is quasi-balanced if and only if it is possible to draw a quasi- 
Eulerian tour going from x 0 to x n and such that, along the lines of the BEST theorem, 
the number of independent quasi-Eulerian tours is given by 


£ = ?i„nhr t -i+'5»,*d ! . (ii) 

X 

where T x is the number of arborescenses, or directed rooted spanning trees pointing at x. 
This number can be expressed as a matrix determinant [36] . To this purpose, we define 
the | AG tv | x \X N \ digraph Laplacian matrix (also called Kirchhoff-Tutte matrix) 


*x,y 



Vx,ye X N 


( 12 ) 


where is the set of states that at time N have been touched at least once by the 
trajectory, that is such that <^4 n + <p£ ut > 0. Notice that Y, x ^ x , y = 0; hence det A = 0. 
Notoriously, the matrix-tree theorem in algebraic graph theory establishes that the 
cofators of the x-th column of the Laplacian matrix give the number of oriented spanning 
trees with root x. In particular, for a digraph we have that 


det A Xn 


(13) 


where A Xn is obtained by removing from A the row and column corresponding to 
Xn. Finally, to obtain the number of distinguishable paths we should identify all 
permutations of edges in the same direction, introducing a Gibbs-type factor. We finally 
obtain 


e|>] 


EW\ 

Fix,!/ Txy- 


det A XN [p] 


n,(^ out -i + 4,x JV )! 

Fix ,y Txy- 


(14) 


where we made explicit the dependencies on the fluxes. This, together with Eq. (ITj) . 
constitutes the main result of this paper. Putting them together, assuming that the 
fluxes ip are large enough so that we can approximate ~ c^° ut = p x and use Stirling’s 
formula ip xy \ ~ exp(p xy \i\ip xy ), we obtain 


P[(p] = det A XN \tp\ e 


(15) 


where 

IM :=!>*„ In-^2- 

x,y T^xyTy 


( 16 ) 
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Interestingly, the digraph Laplacian above can be seen as the generator of a continuous¬ 
time Markov chain for which a particular configuration of fluxes is typical. The same 
prefactor has also been recently derived in Ref. [37] by completely different methods. 


2.3. Large deviations 


It will now be convenient to scale quantities in the extensive parameter N. Let us define 
the flux ratio c p' as the probability that, if a transition occurs, that transition is x <- y, 

Txy _ Txy / 17 \ 

T xy ' at- 1 (f ' ) 

2lx,y Txy N 

Notice that the flux ratio has the meaning of a flux rate, in discrete time. Letting N be 
large but finite (so that the flux ratios remain discrete variables), the probability of the 
flux ratios is given by P'[ip'] = P[Ncp r ] and we define 

mv) 


cm ■■= 




(18) 


where 7:= 7T 


N 

and D( ■ || •) is the Kullback-Leibler divergence. Notice that in 
Eq. (fT5l) while this entropic term decreases exponentially, the determinant increases as 
/VAN but since |An| is at most |X| this term is subdominant f 


■yT'yi 


log P' N W] 


N 


= D(ip' W) + \X, 


N\ 


, hi N In det A XN \ip'] 


N 


N 


(19) 


Then V - - liin/v^ +00 log P' N j iY is the large deviation rate function of the process. We 
thus recover, by a different procedure than usually employed, the Sanov theorem for 
Markov chains [2[ [38; 39] . 

The most probable configuration of fluxes (marked *) is the minimum of the rate 
function, found by solving = 0, which yields cp'* = . This solution is 

indeed unique, as by summing over y we obtain <p' x = Y.y'^xyT’ y i which together with 
T,x T'x ~ 1 grants that Lp' - p is the unique invariant state and ip'* = rz p . 


2-4- Fluctuation relation 

Another testing ground for our formula is the fluctuation relation (see Refs. [10] 40| for 
a general treatment of the fluctuation relation for the currents). Let us reverse the 
direction of all the fluxes by defining 

Ay = <Pyx- (20) 

Notice that in the reversed digraph (A", y>t), x 0 plays the role of the sink and of the 
source. It is obvious that there is a one-to-one correspondence between the number of 
Eulerian tours from x 0 to Xn in (X,(p) and that of Eulerian tours from xn to Xo in 
(A, Then we easily obtain that 

= exp ^ <p~ y In -S3 -. (21) 

x,y Kyx 



§ In the limit N -* oo the flux rates become continuous variables whose probability density scales like 
7V 2 IN, which again provides a subdominant term. 
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Notice that only currents enter this expression. Hence one can finally marginalize to 
obtain a finite-time fluctuation relation for the currents. The functional at exponent in 
the right-hand side is called the entropy production or Lebowitz-Spohn action functional. 

Finally, the fact that there is a one-to-one correspondence between Eulerian tours 
on (X,ip) and of reversed Eulerian tours on (A, yd) implies that the number T Xn of 
oriented spanning trees with root at x n in (X, y?) and that T Xo of spanning trees with 
root at x 0 in ( X , yd) are related by T Xn /T X0 = This is far from obvious as the 

shapes of these spanning trees might be very different. 


2.5. Small fluctuations 


In this section we consider small fluctuations out of the most probable solution, at 
sufficiently large but finite N, evaluating the leading correction coming from the finite¬ 
time term. We set 


<P 


N 



( 22 ) 


where N~ l l 2 5ip « tt p for all edges. Here, it is assumed that “reasonably probable” 
fluctuations scale with %/iV while most probable values scale with N. The leading term 
in the large deviation rate function is of order 0(1). In fact, expanding in Taylor series 
one obtains 




J_ idPxy ^xyTy) 

N Xj y T^xyPy 


E 

x,y 


^Pxy 
xyPy 


(fi^Pxy 'Kxyfo^Py) 


(23) 

(24) 


Let us now show that finite-time corrections are of order N U 2 . Using the log-trace 
formula for the determinant, we have 


det A Xjv [y>] « A |Xjvl det A XN [n p ] exptr log (ld| Xjv | + N 1/2 A XN [n p ] . (25) 

As regards the prefactor, notice that A Xn [tt p ~\ = A^^fTr] diag(p a .) 3 , ;#:XJV . Taking the 
determinant 

det A XJV [7r p ] = det A Xn [-7T] ]“[ p x (26) 

x^xn 

we notice that the determinant of the minor of the Laplacian matrix containing the 
transition probabilities appears. It is well known that, by the matrix-tree theorem for 
Markov chains [41] , we have p XN = A XN \7r]/T(ir) with the spanning-tree polynomial 

rW'EEFW ( 27 > 

T X ttTx 

where the sum is over all possible roots x and over oriented spanning trees T x with 
root x, viz. maximal sets of edges of the graph that contain no cycles and such that 
there exists a unique directed path from each vertex to x. The product is over all edges 
belonging to the spanning tree. 
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Putting all together we obtain 

P[<p\ = N^T(tt) H Px exp {-I[(p] + 6I[ip]} (28) 

X 

where the correction to the rate function is given by 

H[<p\ = N' 1 / 2 tr {A XJV [7T P ] _1 A XJV } . (29) 

Hence the theory allows to obtain a simple expression for the deviations of order iV -1 / 2 to 
the large deviation rate function, when considering a reasonably probable fluctuation of 
size N~ 1 ! 2 . The correction might be simple to compute and contract for the macroscopic 
currents in specific systems such as exclusion processes and the zero-range process. 


3. A ghost field theory for the fluxes 


3.1. Ghosts 


While the dominant behavior at large times is entropic, it is tempting to associate some 
form of energy to finite-time corrections by incorporating the determinant T XN = det A Xn 
into an overall Lagrangian. 

A standard procedure to exponentiate determinants is to introduce anticommuting 
Grassmann variables if - (if x ) x and if = {ifx)x a t ah vertices of the graph. To this 
purpose we notice that det A Xn = det A', where A' is obtained by adding +1 in the row 
and column of A corresponding to xjy, A' xy = A xy + 5 X:XN 5 VtXN , as can be easily seen by 
expanding with Laplace’s formula and keeping into account that det A = 0. Then 

T Xn = det A' = J V(if,if) e^ (30) 


where the integration over the Grassmann variables follows the rules of Berezin’s 
calculus, that we do not discuss here (see Ref. [31] for a review with application to 
matrix-tree theorems). The second expression in Eq. (J5]) is then obtained by identifying 
the LagrangiaiJj] 


£(</>, <P) = ' A [<P] + ifxn i>*n - I[<P\ 


(31) 


At finite times there is a competition between an entropic term and a kinetic term 
of a discretized ghost field of the kind introduced by Faddeev and Popov in Yang-Mills 
quantum theories. 


3.2. Metric and gauge potential 

The analogy to ghost fields can be made stronger. The scope of this section is to recast 
this expression as 

if ■ Aif = (dif, Dif)<$, (32) 

|| An alternative, and possibly more elegant expression, is to write the theory as an expectation value 
m P = f 'D('tp,ip)ip xN ifx N e c ° with respect to the Lagrangian £ 0 (t/>, i/>, ip) = if ■ A[ip]il> - I[tp] 
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with d an appropriate divergence, D = d +A a covariant derivative, A a gauge potential, 
$ a metric. 

First, the proper concept of divergence on graphs is defined on the oriented graph 
G < = ( X, Efl), where we introduce an arbitrary order relation between states x < y which 
establishes an orientation of each edge, that we also denote x < y. We then introduce 
the incidence matrix d with entries 

+1, if y - x 


< 9 , 


y < z,x 


-1, if z = x 
0 otherwise 


(33) 


that fully describes the topology of the graph. Indeed, the incidence matrix is the 
discrete version of the divergence, when acting on scalar fields (defined on the vertex 
set) on the right, and of the gradient, when acting on vector fields (see Ref. j42| for 
elements of discrete differential geometry). 

The latter vector fields are defined as antisymmetric edge observables. Hence 
we move from the fluxes to the activities and currents introduced in Eq. (jHJ), and let 
A' = A + + A - , where Af y S x , y HzTtx ~ T%y Finally we define the positive definite 
metric $ := diag (<Py <z ) y < z and the gauge potential as the \E\ x |A| matrix 


A 


y<z,x 


Tyz/Tyz’ if y = x or 2 = X 
0 otherwise 


(34) 


It is straightforward that A + = <9 T< f><9, where T denotes matrix transposition. Also, 
A~ = d T §A, since 


(d T $A) x , y 


Txy if X * y 

Hz Vyz - Hz Tzy if X = y 


(35) 


We also notice that, introducing the matrix |<9| identical to d but with all -l’s 
turned into +l’s (a sort of unoriented incidence matrix) we have A = A' |<9|, where 
A! = diag This clearly expresses the fact that the covariant derivative is 

not a derivative^! Furthermore, diagonal entries of A' provide an appropriate notion of 
“gauge potential” since, being defined edge-wise and being antisymmetric, they identify 
a discrete vector field. 

Then, we can indeed put the ghost Lagrangian in the form of Eq. (I32|h The metric, 
establishing a vicinity between states, is mediated by the activity, i.e. the total number 
of transitions that occurred between two states. 


3.3. Reflections on the gauge potential 

In Ref. [32] the Author proposed a gauge theory of Markov processes, where the 
underlying symmetry is the invariance of thermodynamic observables under a change 

f From a discrete-algebraic perspective |32] , derivatives are boundary operators. While oriented 
graphs are indeed cellular complexes, i.e. a topological space with proper boundaries and boundary 
relationships such that “boundaries have no boundary”, unoriented ones are not, since their boundaries 
might have boundary. 
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of reference prior, granting the observer-independence of thermodynamics based on 
information theory. There, the gauge potential is a non-fluctuating quantity defined 
by \TL(it X yl'K yx ') modulo gauge transformations. When transition probabilities have a 
physical origin (e.g. they obey the Arrhenius law, the mass-action law, the Kubo-Martin- 
Schwinger condition etc.) the Wilson loops of this gauge potential yield Clausius’s 
expression for the entropy produced along a cyclic transformation. The candidate for 
a fluctuating analogue of this gauge potential would then be ha.((p xy /(p yx ), which is 
obviously different from the entries of matrix A above. However, if we consider small 
fluxes if- y = e~ y we have <p xy = ipf. y + e~ y and 

ln^= ln f 1+ + 2Q L 2 ^ + Q(€ 2 ). (36) 

<Pvx \ <P'xy ' <:nj Txy 

Hence when we consider small deviations from a vanishing configuration of currents we 
do obtain that the vector potential arising from the covariant derivative in the ghost 
field theory is the fluctuating version of the gauge field previously introduced in Ref. [32]. 

It must be noticed that the fact that this identification is approximate might be an 
artifact of the discrete case, and that this discrepancy might fade away for continuous- 
state space diffusion processes. Indeed, the large deviation rate function for diffusion 
processes with empirical measure /i and current j is given by mm 

= \ f U -j -jn)dx (37) 

on the assumption Vj = 0, where j y = Ffi-DV/i, D being the diffusion matrix and F the 
drift vector. Notice that in this expression again a fluctuating symmetric property /j D , 
with the meaning of dynamical activity, plays the role of the metric. The gauge theory 
for diffusion processes was analyzed in Ref. [45], an d there indeed the vector potential 
is locally linear in the currents as in Eq. (13611 . and not a logarithmic expression; its 
candidate fluctuating counterpart would then be A - (/ilf) _1 j. Therefore, generalization 
of our results to diffusions in continuous space, e.g. along the lines initiated in Ref. [46]. 
might indeed lead to an elegant and exact gauge theory for finite-time fluctuations. 
A rigorous treatment of the small-step limit from discrete to continuous diffusion is 
postponed to future inquiry. Techniques based on fermionic fields, similar to the ones 
employed in this paper, for Fokker-Planck equations have been employed in Ref. [43] . 

4. General Markov jump processes 

In this section we generalize the results of Sec. [2] to Markov jump processes. As we will 
see, application of the BEST theorem leads to no significant variation on the theme, 
while some care has to be paid in handling the waiting times spent at states, which 
enter the game in a crucial way. Below we consider two strategies to deal with them: 
moving to the Laplace representation (as in Refs. [151. 1 47]. among others), which is quite 
straightforward but requires an awkward (still, feasible) inversion, or considering the 
joint probability of the fluxes and of the so-called empirical measure. Again, we re¬ 
obtain known results from large deviation theory, putting some emphasis on certain 
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subtle aspects related to the statistics of total number of jumps performed up to some 
time, that we feel have been overlooked. 


4-1. Jump processes in Laplace space 


Let U(t) = exp[f(bb - 11)] be the propagator of a stationary ergodic continuous-time 
Markov semigroup on a finite state space. Entries of the off-diagonal part W of the 
generator are the positive rates w xy at which transitions x <- y occur, with units of 
inverse time. Entries co y = Y,x w xy of the diagonal part 11 are exit rates out of states. 
We will work in Laplace representation, assuming it exists (the connection between 
semigroups, the resolvent formalism and the Laplace transform is ruled by the Hille- 
Yosida theorem, see [38] for a simple introduction). For all u > 0, 

U(u) = J e~ UoI U{t)dt = (ul + H - IE) -1 , (38) 

where I is the identity. It can be proven that ||IE|| < ||u;/ + H|| with respect to the usual 
operator norm, under which condition the following operator geometric series can be 
produced 

Too 

(H - W)- 1 - ll- 1 X! (LLH- 1 )^. (39) 

iV=0 

Making the matrix products explicit, Laplace transition amplitudes can be expressed as 
sums over paths 


N x 


with path weight 

P N [x](u) 


n 




(jJ + UJn i<N UJ + UJi 


(40) 


(41) 


where path x = (xi) performs N jumps between the constrained end-points xq and 
xn = x t . For notational simplicity, wherever possible we replace state X; by i. 

The analysis then follows as for restless Markov chains. We introduce the fluxes <p 
as in Eq. (J2J) to obtain 


P N \v]{u) 



detA ^] e -W]M 

U) + 0J]\f 


(42) 


where 


I N [<p](uj) := X! Txy hi 

x,y 


Wxy^Py 


(43) 


An important note is here in order. The probability of the fluxes Eq. (j42T) is labelled by 
the total number N of jumps to account for the fact that the possible values that these 
fluxes can take depends on N, since Y, x ,y Txy = N. In mathematical terms, the filtration 
of the process (the time sequence of the possible events) is itself a stochastic process 
regulated by the probability of N jumps up to time t. We will return on this issue later. 
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While Eq. (14 2 j) is in Laplace parameter, it can already be employed in practical 
experiments. In fact, if instead of halting the Markov process at a given time t we 
conduct experiments in variable time sampled with probability P(t ) = uie^, then the 
probability distribution of the fluxes for such protocol is given by Pn,w [<P] = uP N [(fi](u). 
In this case, since typical trajectories halt at time a; _1 and prolonging this time makes 
the probability smaller, large fluctuations are very rare and it is necessary to retain the 
full structure of the probability density function. 


4.2. All-time probability of the fluxes 

If instead we insist on performing experiments in a given time, we need to perform the 
inverse Laplace transform. Isolating the explicit dependency on u we have 

1 


p N[ip](t) = det A Xt [gp] e Jj vbd(°)£ 1 


n 


( 1 +uj/cu x y* 

The inverse Laplace transform to be calculated is given by 

1 


(t). 


C 


-1 


n 


(t) = Y l Qx( t ) ( 


I Mx t 


(44) 


(45) 


(u x + oj)** 

where Q x {t ) is a polynomial of degree <p x -1 in t. The general idea is that Y\ x {u x + oj)^ x 
is the common denominator for denominators (uj x + u) mx , for all n x = 0,... ,(p x - 1. If 
there exist coefficients C x x such that 

1 „ , 1 


n 


= E E or 1 - 


(u x +u)V* x 1 * (uj x +u) mx 
given that the Laplace transform of these latter terms is known, 


C 


-1 


(uj x + U)) r 


(t) = 


f-TLx 1 


-w x t 


(n x - 1)! 


(46) 


(47) 


we obtain 


Q*(t) = E c ?— 

m x =0 


Ux 


The task is then to find the coefficients C™ 
Appendix B, 


C nx = (_ 1 )V>x-l-n a £ Yl 


(48) 

We derive the following formula in 


^ ^q (Tv/q)' 


{m y } y*x (u y - U x )Vy+™ y { ^ } I = I mq \ 
Uly+x m y 

— ( -Px~Tl x ~^- 


Eji“ 1 om q 


(49) 


Though not very practical, this formula is programable in software simulations. 


4-3. Joint fluxes and empirical measure 

The empirical measure accounts for the resting time at states. In this section we derive 
an exact formula for the joint statistics of the fluxes and of the empirical measure at 
time t, for a given number of jumps N. 
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Before jump Xi -* x i+ i, the process waits an interval of time t t at state x i} for a 
total elapsed time b; = t. We then introduce the empirical measure as the fraction of 
time spent at x: 

r x [x,t] := Yjtibhx- ( 50 ) 

i 

Notice that Y,x T x[ x ,t] = £,b = b We can then express Eq. (ITTj) as 

P N [x](co)= f Vt e~ iuJ+UJN)TN Y[w i+hi e-( Ui+oj)ti (51) 

J i<N 

= dt e~ wt Vt P N [x,t](t) (52) 

where all integrals range from 0 to +oo and 

P N [x,t](t) := 6(t- ^bV^^ n e ~ WiU ( 53 ) 

' i 7 i<N 

= 5(t-Y j r x [xfl]y-^ u ^ t ^Ylwf; v[x] . (54) 

x x,y 

is the probability that up to time t, given N jumps, trajectory x,t has been performed. 
We now want to marginalize for <p,T. To marginalize out the waiting times, noticing 
that the above expression only depends on t[£E,£], we need to calculate 

[vt = (ss) 

J X X H^x- 

Notice that for large enough ip the factorial p x \ is counterbalanced by that coming from 
the Eulerian tour counting. We then obtain, using Stirling, 

PN[<P,T](t) =<$(*-£ Tx) det A Xn exp r] (56) 

X 

where 

InW, r] := Yj Txy In ^ xy +Y,V v x - (57) 

x,y 'Wxy'Ty x 


f.f. Large deviations of the fluxes 

When considering large deviations one requires that the fluxes and the empirical measure 
grow linearly in time 

<p ~ t(p (58a) 

r ~ tr, (58b) 

where the dot serves as a notation. The joint rate function for the fluxes and the 
empirical measure has been discussed in Ref. |2£) and, interestingly, it can be further 
contracted for the joint currents and empirical measure |44j[49j. As a small improvement 
we can work out an expression for the large deviation rate function of the fluxes. We 
define 


7jv(0 t 1 I N (t-). 


(59) 
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Notice that we retain the explicit dependence on N for reasons that will be explained in 
the following section. Let us take the variation of In with respect to r, accounting 
for the normalization E x f x - 1 by means of a Lagrange multiplier A. Setting 
8/5f x [/jv + - 1)] = 0 yields f* = fl x /(ui x + A), where the multiplier satisfies 


E 


<Px 


uj x + A[0] 


- 1. 


We then obtain 


In[<p] = E Txy ln 


x,y 


ip X y(UJy + A[y>]) 
^xy^Py 


E 




u x + A [<p] 


(60) 


(61) 


This, together with the constraint Eq. (16011 provides an expression for the large deviation 
rate function for the fluxes, with an important caveat that we next discuss. 


4.5. Total number of jumps and large deviations 

Derivations of the large deviation rate function for the fluxes and the empirical measure 
often employ the smart trick of comparing the probability of a fluctuation with that 
of a system where that configuration is typical. The rate function found by the above 
mentioned authors is (see Eq. (11.7) in Ref. (44]): 

J((p,f)= X! Txy In -S-SE + Y j t x u x ~Y, Tx,y (62) 

x, y ^xy^y x x,y 

= i N (<f,T)-j (63) 

The second line identifies J with our expression but for an additional term -N/t. From 
our perspective, the origin of this term is subtle and we haven’t seen it treated elsewhere. 
In the following we attempt a discussion. 

First we notice that there is no such issue in the case of restless Markov chains, 
because N - t with certainty and additive constants to the rate function are irrelevant. 
For jump processes, the total number of jumps At is a stochastic variable itself and it 
obeys a large deviation principle. It might be that at some time t a number N oc t 
of jumps have occurred, but (since we are interested in large deviations!) we need 
to keep into account that also N might fluctuate. In other words, the two extensive 
parameters N and t are decoupled and one needs to specify how one scales with the 
other. Importantly, while <f> determines N, its statistics depends on N in a subtle way. 
In fact, a Markov jump process is built by superimposing Poisson processes on top of a 
Markov chain [50], so that the sequence of measurable values that the fluxes can take 
(the so-called filtration of the process) is itself a stochastic process. 

The simpler way out of this riddle is to just suppose that N = co*t independently 
attains the minimum of its large deviation function, that is, that we look at the most 
probable fluxes assuming that the total number of jumps up to some time is the most 
probable one. Then in Eq. (163|) we obtain an additional constant that does not affect 
the large deviations. A posteriori , it is intuitive that the most probable configuration 
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of fluxes should be such that, of all the jumps out of y, the ratio of those that go to x 
should equal to the probability that the next jump out of y is to x: 


V*y 


W 


xy 


UJr,, 


(64) 


Notice that w xv /u y are transition probabilities of a restless Markov chain. Then 


J xyl y 

fl* = to* p x , p being the steady state of the restless Markov chain 

check that this is the stationary solution corresponds to A = 0 and 

/ \ -i 

Px 


CO 


= £- 


to x 


One can easily 


(65) 


As a consistency check, notice that these are indeed the exact expressions when one 
takes all identical exit rates co x = to*. Therefore, we do obtain the usual rate function 
as in Refs. mmm, when the total number of jumps is assumed to attain its most 
probable value. However, the assumption is quite a strong one, as it has been shown that 
there exist kinetically constrained models where the total number of jumps undergoes a 
dynamical phase transition between an active phase with N = to*t and an inactive phase 
with N subextensive in time [5T| 152) . 

We hold that a more complete and correct treatment should instead keep into 
account that the total number of jumps fluctuates, and that this might imply 
modifications of the formula usually found in the literature. The large deviation rate 
function of the total number of jumps has been the subject of a prolific branch of studies 
[7711 [52] . For the case of all identical exit rates to x = to*, the probability of N total jumps 
is Poissonian with average to*t 

P[iV](t) = e-" ,, h^A (66) 

with rate function 

■ N 

I[N]=co* + N\n—. (67) 

to* 

This piece of information should somehow feed back into the rate function of the fluxes 
if the general question to be asked is what is the probability of a configuration of 

fluxes, regardless of the total number of jumps. We demand this task to future inquiry. 

Departure from Poissonian statistics in a low-density fluid, treated with linearized 
Boltzmann equation, is found in Ref. [53]. Analogies and differences between large N 
and large t ensembles have been studied in Ref. [53]- 


5. Conclusions 

In this paper we established a fascinating connection between Large Deviation Theory, 
combinatorics, and certain techniques and concepts from the Quantum Field Theory of 
Yang-Mills theories. The main result is a generalization of known results [29j 01 , 09 ] 
regarding the large deviation rate function of the fluxes and of the empirical measure 
for Markov processes (level 2.5 large deviations) to the full probability distribution. 
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The result relies on a slight generalization of the so-called BEST theorem in graph 
combinatorics to Eulerian tours with open ends. On the other hand, since Eulerian 
tours can be counted by spanning trees, recently proposed methods to express matrix- 
tree determinants as Fermionic field theories are applicable. The Lagrangian so found 
displays an interaction of the ghost field with an underlying metric, representing the 
activity of the system, and with a gauge potential that is related to the nonequilibrium 
gauge potential that the Author explored in a previous publication [35]. 

The mathematics of spanning trees is very rich and full of connections, e.g. to knot 
theory via the Alexander’s polynomial [55], to spin networks in Loop Quantum Gravity, 
to Feynman diagrammatics. A simple open question that already emerges is whether our 
slight generalization of the BEST theorem to quasi-balanced digraphs might generalize 
to any digraph, by a suitable decomposition in terms of several Eulerian tours with 
open ends. Finally, establishing the statistics of the currents at all times might play an 
important role in the study of finite-time thermodynamics. 

An open problem that will deserve further attention is that of the apparent 
dependence of the large deviation function, as we derived it, on a given total number of 
jumps and, more in general, a more careful characterization of the statistics of the total 
number of jumps, which is affected by the stochastic trajectory and in turn determines 
the set of measurable events. 
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Appendix A. BEST theorem for quasi-balanced digraphs 

Let ( X , q?) be a connected balanced digraph with vertex set X. Let ip xy be the number 
of arrows y -> x and ip£ ut = Y, y Pyx be the out-degree of vertex x. 

Definition 1. A quasibalanced digraph is a digraph such that one vertex x init has one 
more outgoing arrow than ingoing ones, one vertex x&a has one more incoming arrow 
than outgoing ones, and all other vertices are balanced. 

Definition 2. An open Eulerian tour on a quasibalanced digraph is a path that starts 
at Xi n i t , ends at x& n and traverses each arrow exactly once. 

Theorem 1 . A connected quasibalanced digraph admits an open Eulerian tour. 

Proof. By adding an arrow Xi n i t x& n one obtains a connected balanced digraph, which 
admits a closed Eulerian tour starting and ending at £i nit . By cyclically permuting the 
arrows in this tour one can obtain a Eulerian tour whose last arrow is Xi n it ■«- Xfi n . 
By removing this arrow one obtains an open Eulerian tour on the quasibalanced 
digraph. □ 
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Theorem 2. The number of open Eulerian tours on a connected quasibalanced digraph 
is given by 

(A.l) 

X 

where T XRn is the number of oriented spanning trees with root in x fin . 

Proof. We follow the arguments found in Stanley’s book [35], pp. 54-57, for balanced 
digraphs, with slight modifications. Consider an open Eulerian tour, that by the 
previous theorem exists. For all vertices x + xg n let e(x) be the last arrow that exits 
x. Then the set T = { e(x),x * Xg n } is an oriented spanning tree with root in x fin (see 
[35] . Proof of Claim 1). Conversely, given a spanning tree with root in xg n one can 
construct open Eulerian tours starting at x init and picking random consecutive edges 
that do not belong to the spanning tree, unless there are no edges left out, in which 
case one will pick an edge e(x) e T and exit x for the last time. We need to show that 
by this procedure one doesn’t get stuck at Xg n before traversing all other arrows, that 
is, that we do not need to “lift the pencil”. If that was the case, then it means that at 
least a simple cycle of arrows C (with no crossings or multiple edges) has been left out 
of the tour 7 so far generated (since E \ 7 is a balanced digraph). Since the spanning 
tree T does not contain cycles, one of these arrows e € C out of a vertex y is not in the 
spanning tree, and e(y) does not belong to the cycle. Then either we exited y by e(y) 
before using e, which is in contradiction with the above procedure, or e(y) has not been 
used, which implies that it belongs to another simple cycle C' that has been left out of 
7 . Along this other simple cycle there exists an edge e' i T, and one can continue the 
reasoning to obtain a contradiction. Then there are fix (Ax Ut - 1 + 5 X)Xfln )! ways to pick 
such random arrows, with the -1 accounting for the edges in the spanning tree, but 
for vertex xg n . Since every choice of a tree T and every random selection of remaining 
arrows generates a different Eulerian tour, and every Eulerian tour can be so described, 
we can conclude. □ 


Appendix B. An inverse Laplace transform 


The task is to find the inverse Laplace transform of 

1 


F(cu) = n 


(B.l) 


(oj+u x y* 

The result is a generalization of Heaviside’s expansion theorems [56] • We employ the 
well-known residue formula 

fit) := /r 1 [F(u)] (t) - Y , Res -U) X ]. (B. 2 ) 


Since -u x is a pole of degree ip x , we have 


m =e 


lim G 


ifPx 1 ) 


(<Px - 1)! 


(a) 


(B.3) 
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where we introduced 


-j^[ eUlt EK^ + UJ v) 


Vy 


y±x 


(B.4) 


Interpreting as a moment generating function, instead of directly evaluating 

moments it’s simpler to evaluate the cumulants 

r(<?)/ 


■- —In G^(uj) 

d<jj q 

= tS, A + (-!)«(? - 1 )!£ 


Ty 


y±x ^y 


t8q t 1 + Edx q \(v). 


(B.5) 

(B. 6 ) 

(B.7) 


The generalization of the chain rule to higher derivatives is given by Faa di Bruno’s 
formula, that simplifies for the exponential function g4 0 ^(w) = exp H x °\cu): 

,, n\ p r • " ma 

ci-vo=y m /• „, n 




m\ \... m p \ [ q\ 

= G ( °\u)B p u), ..., Hi P \co)) 

where the sum runs over all rri \,..., m p > 0 such that 

p 

E 2 m 9 = P- 

9=1 

The second line identifies the complete Bell polynomials. Let us now isolate powers of t 
using of the binomial formula, and rearrange summations. After some work we obtain 

g?>( w ) = g‘°V) y [t +j?< 1) ( w )]"" - / ! m , n i H * (w) 


(B. 8 ) 

(B.9) 

(B.10) 


mfl.... m p \ q = \ [ q\ 


= y y 


mi 


m\\ 


Wnt^o ra^mr-n*)! 




p! 


rUil-.-mp! P L 






= gS°V) yyy 




pi 


n 


H ( x q \u) 

q\ 


where we have defined fhi = mi - n x , ffijy i = and now the sum over the {rhj} is 
constrained by 

p 

E Qkhq =p- n x . (B. 12 ) 


(B.ll) 


Finally 


9=1 


V^-l fn x 

m = y«-"-‘ y — c? 

i n x =0 ^x- 


(B.13) 




y*x 


- 1 1 ^ 

{^} 9=1 m 9 ! Uy*x K-W,) 9 


(B.14) 


with 
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We can further introduce a multinomial expansion with Y, y *x m q,y ~ m q, define m y = 
Y, q qm q ,yi replace £{m 9i „} by T,{ my } keeping into account the constraints 

Hy*x m y = (p x -n x -l and Y,q= i _1 qm q , y = my, to obtain 

(‘ Py/q) mq ’ y 


c%* = (-i)—n( 


U),. - u. 




y±x 


E 


^x~y i 

\-vv y y n n_ 

(fei {m^y} 5=1 y*x m K - ^) qmg ’ : 

1 E nn^EXfB.is) 


{m v} n y*x(uy - u x )^y +m y { flfl q y} y “ "W 


Appendix C. Example 


Consider as the state space the following graph: 


a - b 



A configuration of fluxes identifies the following quasi-balanced digraph: 


(C.l) 



The Laplacian matrix reads 


/ 

2 

0 

0 

-1 

\ 


-2 

3 

0 

-1 



0 

-1 

1 

0 


V 

0 

-2 

-1 

2 

/ 


(C.2) 


(C.3) 


Let us enumerate the 6 independent Eulerian paths going from the source a to the sink 
d (to do so systematically one can proceed as follows: start at a, and every time choices 
ramify into several possible directions, choose the one with the lowest letter until all 
transitions are exhausted; then choose the second lowest at the latest ramification, and 
so on...): 


abcdabdbd 

abcdbdabd 

abdabcdbd 

abdabdbcd 

abdbcdabd 

abdbdabcd (C.4) 

Hence, accounting for the Boltzmann factorials for the fluxes, we obtain 6 x 2! x 2! = 24 
Eulerian tours. 
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There are two spanning tree pointing at d: 


a ->- b a -»■ b 

(C.5) 

d ■* - c d -«-c 

The first occurs 2 times because of the choice of edge a -*■ b, the second 4 times because 
of the choice of edges a -*■ b and b -*■ d, yielding 6 spanning trees (that this number 
coincides with the number of independent paths is a mere coincidence). Finally we have 
that this number has to be multiplied by the factorials at the vertices, (3-1)! at b and 
2! at d, yielding 24 Eulerian tours. 

The cofactor reads 



det A d = 


2 0 
-2 3 

0 -1 



(C.6) 


and indeed it yields the number of oriented spanning tree rooted at d. 
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